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Abstract Energetic nonthermal particles (cosmic rays, CRs) are accelerated 
in supernova remnants, relativistic jets and other astrophysical objects. The 
CR energy density is typically comparable with that of the thermal compo- 
nents and magnetic fields. In this review we discuss mechanisms of magnetic 
pl^ I field amplification due to instabilities induced by CRs. We derive CR kinetic 

and magnetohydrodynamic equations that govern cosmic plasma systems 
comprising the thermal background plasma, comic rays and fluctuating mag- 
netic fields to study CR-driven instabilities. Both resonant and non-resonant 
instabilities are reviewed, including the Bell short- wavelength instability, and 

O , the firehose instability. Special attention is paid to the longwavelength insta- 

bilities driven by the CR current and pressure gradient. The helicity produc- 

^ I tion by the CR current-driven instabilities is discussed in connection with 

the dynamo mechanisms of cosmic magnetic field amplification. 
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1 Introduction 

Acceleration of cosmic rays (CRs) in the Galaxy by the first order Fermi 
mechanism is believed to be very efficient. Most of the theoretical studies of 
shock acceleration agree on its potential to convert, under favorable condi- 
tions, 50% or more of shock mechanical energy into the CR energy. Observa- 
tional estimates of the supernova remnant (SNR) shock power require, on the 
average, a 15-30% conversion efficiency to mainta in the observed CR energy 
again st losses from the Galaxy (see, e.g.. iBcrczinskii et al.lll990llDrurv et al.l 
Il989[ ). However, this acceleration mechanism is fast enough only if it is self- 
sustained; accelerated particles must be scattered across the shock at an 
enhanced rate (to gain energy rapidly) by magnetic irregularities amplified 
by the particles themselves. Relying on the background magnetic irregulari- 
ties (interstellar medium [ISM] turbulence) would result only in a very slow 
acceleration. 

Fortunately, freshly accelerated CRs indeed comprise enough free energy 
to d rive plasma ins tabilities thus bootstrapping their own acceleration (see, 
e.g., IZweibellll979[ ). While they are accumulated in a relatively thin layer 
near a shock front, their pressure gradient is built up. Furthermore, they 
stream through the inflowing plasma so that their pitch-angle distribution is 
anisotropic. They also provide an electric current and induce a return current 
in the upstream plasma. 

Instabilities driven by the above sources of free energy may loosely be cat- 
egorized as follows. First, an ion-cyclotron type, resonant instability (driven 
by the CR anisotropy) amplifies Alfven and magnetosonic waves, with no 
major changes to their dispersive properties and the macroscopic state of the 
medium near the shock. However, the amplified waves make the CR pressure 
and current to build-up rapidly through an enhanced CR scattering and en- 
ergy gain. Second, there is a non-resonant firehose type instability driven by 
the CR pressure anisotropy. In contrast to the resonant instability, the fire- 
hose instability changes the Alfven wave dispersive properties by making the 
growing mode aperiodic. So does the current driven non-resonant instability. 
The renewed interest to this instability has been sparked bv I Belli ( 20041 ). who 



revealed its potential to strongly amplify the background magnetic field. In- 
deed, a formal analytic solution in which the instability driver is balanced by 
the nonlinearity indicates t hat the ins tability saturates only at very h igh ani- 
plitudes, SB > Bp (see e.g.lBeU and L ucck 2001, Bell 200 5.lMarcowi th et alj 
120061 . ICaprioh et al.ll2008l IVladimirov et al.ll2009l iMalkov et al.Tl2oi2bD . Fi- 



nally, the CR pressure gradient in the shock precursor drives acoustic pertur- 
bations. All these instabilities should be treated on a unified basis, as they 
arc driven by the anisotropic inhomogeneous CR plasma component near a 
shock front. An attempt of such treatment is presented below. However a 
complete nonlinear study of these phenomena is a formidable task, yet to be 
accomplished. 

While the above instabilities, clearly associated with coUisionlcss shocks, 
will be central to the present review, CRs are also known to drive instabilities 
crucial to their confinement regardless of the way they are accelerated. For 
example, a sufficiently dense CR cloud released into the ISM will drive Alfven 



wav es which, in turn, w i ll scatter the CRs , thus delaying their e scape (see 
e.g. iPtuskin et alJ l2008l lOhira et all l201ll iMalkov et all l2012al lYan et al.l 
|2012[ ). Moving further out to the CR confinement in the galaxy, the so-called 
Parker instability is known to be important, in addition to the Alfven wave 
self-generation by escaping CRs. 

The diffusive shock acceleration (DSA) mechanism is based o n repeated 
shock crossings with a ~ Us/c particle energy gain per cycle (see iKrvmskiil 
1977, 'BcU"l978', 'Bland ford and Eichled [TOSTJ iBerezhko and Krvmskiall988l 
Jones and Ellison 1991|) . While doing so, particles diffusively escape from the 



shock up to a distance Lp ^ n (p) /ug. Here k is the momentum dependent 
diffusion coefficient and Ms is the shock velocity. One should expect then 
an extended (~ Lp) shock precursor populated by accelerated protons and 
electrons so that synchrotron radiating electrons may make it visible. High- 
resolution X-ray observations have revealed thin X-ray synchrotron filaments 
and fast evolving clumps in synchrotron emitting supernova shells. The fila- 
ments arc much thinner than Lp because the TeV regime electrons are con- 
fined in a narrow layer around the shock. Most likely they are limited by fast 
synchrotron cooling due to the X-ray emission in a highly amplified magnetic 
field (see for revie w ICassam-Chenai et al.ll2007l lRevnoldsll2008llVinkll2012l 
iHelder et al.ll2012l ). The synchrotron emission clump s with a year time scale 
variability observed with Chandra observatory by lUchivama et al.l (|2007| ) 
can be associated w ith strong intermittency of the amplified magnetic fields 
(JBvkov et al.ll2008l ). Moreover, a quasi- regular set of strips of synchro tron 



emission resolved with Chandra in Tycho's SNR bv lEriksen et al.l ( 201lh po- 



tentially can be used to study a specific angular dependence and the spectral 
properties of nonli near mechanisms of magnetic field amplification by CR- 
driven instabilities ( Bvkov et al.ll2011al) . 



According to the widely accepted view, the particle diffusion coefhcient 
K should be close to the Bohm value, k ~ erg (p) /3, which requires strong 
magnetic fluctuations 6Bk ~ Bq at the resonant scale k ^ \/rg{p). The 
high level of fluctuations is achieved through one of the instabilities driven 
by accelerated particles. A number of CR driven instabilities have been sug- 
gested to generate magnetic field fluctuations. The first one is the well known 
ion cyclotron resona nt instability of a slightly anisotropic (in pitch angle) CR 
distributio n (sec c.g. lSagdeev and Shafranovlll96lllZweibellll979LISchlickeiseil 
I2OO2I [Amato 2OII1). The free energy source of this instability is potentially 
sufficient to generate magnetic field fluctuations needed to s catter CRs ahead 
of the shock (see e.g. lBel]|[l97llMcKenzie and Voelklll982[ ). 

{SB/B^f ^ MaP'^/puI (1) 

where Ma ^ 1 is the Alfvenic Mach number, P'^^ is the CR pressure, p 
is the gas density and Ug is the shock velocity. However, the actual tur- 
bulence level was shown to remain moderate, SB ~ Bq as this is a reso- 
nant kinetic instability that is usually suppresse d by a quasilinear isotropi- 
sation or particle trapping effect s easily (see e .g. iMcKenzie and Voelklll982l 
lAchterberg and Blandfo^ll986l . IZweibe]|l200l . 

The second instability, is a nonresonant instability driven by the CR cur- 
rent. The advantage of this instability seems to be twofold. First, it cannot 



be stabilized by the quasilinear deformation of the CR distribution function 
since in the upstream plasma frame the driving CR current persists, once the 
CR cloud is at rest in the shock frame. Second, it generates a broad spectrum 
of waves, and the longest ones were claimed to be stabilized only at the level 
SB ^ Bn. due to the lack of efficient stabilization mechanism at such scales 
(see e.g. lBel]||2004). Within the co ntext of t he CR acceleration, t his instability 



was studied by lAchterbera (|l983l) (see also IShapiro et alj 19981). but the fast 



regi me of the nonresonant instability was found by iBell and Lucekl (|200l[ ) 
and lBelll (l2004f). and therefore the instability is often referred to as Bell's in- 
stability. [Bell (J2004f ) pointed out that in the instability is driven by a fixed CR 
return current through the Ampere force jcr x B. It should be noted, however, 
that the dissipation of the return current due to the anomalous resistivity 
still needs to be addressed . The effect of a finite plasm a temperature on the 
instability was studied by IZweibel and EverettI ( 2010f ) . Actually, as we will 



show below, both the resonant and the Bell instabilities are interconnected, 
they are driven by the CR drift relative the background plasma. Moreover, 
in the case of the modes propagating along the mean magnetic fields the two 
instabilities are simultaneously influencing the same modes. The dispersion 
relations of the modes are strongly influenced by the presence of the CR cur- 
rent are markedly different from the standard MHD modes. The dispersion 
relations of the modes strongly influenced by the presence of the CR current 
arc markedly different from the standard MHD modes. The dispersion rela- 
tion in the longwavelength regime (where the mode wavelengths arc larger 
than the bulk CR gyroradii) can be also strongly modified by the pondero- 
niotive forces in duced by Bell's turbulenc e. The longwavelength instability 
has two regimes ( Bvkov et al.ll2011a 120121) . The first regime is prominent in 



the intermediate range where the mode wavelength is above the CR gyroradii 
but below the CR mean free path. It is discussed in ^4.41 and is associated 
with a dynamo type instability driven by the nonzero helicity, which is, in 
turn, produced by the short scale CR-driven turbulence. The intermediate 
wavenumber range is rather narrow in the case of the Bohm-type CR dif- 
fusion. The modes with wavelengths larger than the CR mean free path 
are subject of non-resonant long-wavelength instability caused by the pon- 
deromotive force acting on the background plasma that is induced by Bell's 
turbulence. We discuss the long wavelength instability below in 



The third instability is an acoustic instability (also known as Drury's 

i nstability) driven by the pressure gradient of accelerated CRs upstream 

(iDorfi and Drurvlll985l . [Drurv and Falldl986l[D^urv and DowneJ2012l[Schure et al.l 
|2012[ ) . The pressure gradient is clearly a viable source of free energy for the 
instability. So, among the macroscopic quantities varying across a strong 
shock, the pressure jump is the most pronounced one in that it does not 
saturate with the Mach number, unlike the density or velocity jumps. 

The acoustic instability has received somewhat less attention than the 
first two. Moreover, in many numerical studies of the CR shock acceleration, 
special care is taken to suppress it. The suppression is achieved by using 
the fact that a change of stability occurs at that point in the flow where 
dhiK/dhip ~ — 1 (for both stable and unstable wave propagation direc- 
tions, of course, if such point exists at all). Here p is the gas density. Namely, 



one requires this condition to hold identically all across the shock precursor, 
i.e., where the CR pressure gradient VP*^'' 7^ 0. Not only is this requirement 
difhcult to justify physically, but, more importantly, an artificial suppres- 
sion of the instability eliminates its genuine macroscopic and microscopic 
consequences, as briefly discussed below. 

Among the macroscopic consequences an important one is the vortic- 
ity generation through the baroclinic effe ct (missalignment of the density 
and p ressure gradients Vp x VP ^ 0, e.g. iRvu et all (jl993l ). iKulsrud et al.l 



( 19971 )). Here VP may be associated with a quasi-constant macroscopic CR- 



gas pressure gradient VP^'', generally directed along the shock normal. Vari- 
ations of Vp are locally decoupled from P'^'', unlike in the situation in a gas 
with a conventional equation of state where P = P (p) and where the baro- 
clinic term vanishes. The vorticity generation obviously results (just through 
the frozen in condition) in magnetic field generation, so that the field can be 
amplified by the CR pressure gradient. More importantly, this process ampli- 
fies the large scale field, required for acceleration of high energy particles. Fur- 
thermore, the amplification takes place well ahead of the gaseous subshock. 
The both requirements are crucial for improving high energy particle confine- 
ment and making the shock precursor shorter, in agreement with the observa- 
tions. Large scales should be present in the ambient plasma as a seed for their 
amplification by the acoustic instability and could be driven (or seeded) by 
wave packet modulations. Apart from that, they result from the coalescence 
of shocks generated by the instability, and fro m the scattering of Alfven wave s 
in fc-space by these shocks t o larger scales Malkoy and Diamond! (J2006[) , 



[Diamond and Malkovl (|2007t) . iMalkov and Diamond! (|2009t) . Note that the 
Bell instability is essentially a short scale instability (the maximum growth 
rate is at scales smaller than the gyro-radii of accelerated particles). At larger 
scales the magnetic field growth rate is dominated b y the modified reso - 
nant and the longwavelength nonresonant instabilities (JBvkov et al.ll2011a l. 
ft should be noted that vorticity (and thus magnet ic field) can be efficiently 
generated also at the subshock fsee e.g. | McKenzie a nd Westphal 1970, Bykoy 
1982lll988llKevlahan|ll997llKulsrud et al.lll997llGiacalone and Jokipiill2007f 



Beresnvak et al.ll2009 



Fraschettill20 1 31 ) . This would be too late for improving 



particle confinement and reducing the scale of the shock precursor. A more 
favorable for acceleration scenario is the above discussed field amplification 
in the CR shock precursor. 

Now the question is which instability dominates the CR dynamics? Given 
the finite precursor crossing time, it is reasonable to choose the fastest grow- 
ing mode and consider the development of a slower one under conditions 
created by the fast mode after its saturation. The Bell instability is likely to 
be efficient at the outskirt of the shock precursor where the CR current is 
dominated by the escaping CRs of the highest energies. The pressure gradi- 
ent and the pitch angle anisotropy are strong enough to drive the acousti c 
and resonant instability in the shock precursor (see e.g. lPelletier et al]l2006() . 
Recall that the anisotropy is typically inversely proportional to the local 
turbulence level which is usually decrease with the distance from the shock 

Within the main part of the shock precursor, both the CR-pressure gradi- 
ent and CR current are strong, so that the nonresonant CR-driven instabili- 



ties are likely to be the strongest candidates to govern the shock structure. In 
fact, these instabilities are coupled, not only by the common energy source 
but also dynamically. But first, it is important to identify conditions under 
which one of the instabilities dominates. 



2 Cosmic plasmas with cosmic rays: the governing equations 

In this section we discuss the governing equations for MHD-typc flows of 
a cold background plasma interacting with cosmic rays. In most cases the 
cosmic ray particles are not subject to binary Coulomb or nuclear inter- 
actions with the background plasma particles. The interaction between the 
two components is due to both regular and fluctuating electromagnetic fields 
produced by the CRs. The momentum equation for the background plasma, 
including the Lorentz force associated with these fields is given by 

/ ^ — \ 1 

pi-^ + (uV)a j = -Vpg + -J X B + e (np - Ue) E, (2) 

where B is the magnetic field induction, E - the electric field, u - the bulk 
plasma velocity, pg - the plasma pressure, j - the electric current carried by 
the background plasma. We assume quasi-neutrality for the whole system 
consisting of background plasma protons of number density rip, electrons of 
number density n^, and cosmic rays of number density ricr- For simplicity 
we consider cosmic-ray protons only such that rip + ricr = ^e, and typically 

'lev ^^ ' ^p ■ 

The magnetic field is assumed to be frozen into the background plasma 



E=-i 
c 



ux B 



(3) 



Both the background electric current j and the electric current of acceler- 
ated particles j'^'" are the sources of magnetic fields in Maxwell's equations, 
where the Faraday displacement current was omitted for the slow MHD-type 
processes 

V X B - ^ (j + r) . (4) 

Then, for the quasi- neutral background plasmas, using Eq.Q, Eq.([3]) and 
Eq.Q, one can write the induction equati on and the equation of motion o f 
the background plasma in the form used bv lBelll (|2004 ). lBvkov et aP (|201ld ). 
ISchureandBeiil(|201lD 



3T\ ■— 

— = V X (S X B), (5) 

p ( ^ + (GV)a') = -Vpg + i-(V X B) X B - i(r - en,,S) x B. (6) 



The microscopic CR-dynamics can be described by a kinetic equation for 
the single-particle distribution function / that has the form 

where t he CR particle energy is f , is the momentum rotation operator 
(see e.g. lToDtvginlll983llBvkov et al.ll2012[ ). There are no Coulomb collisions 
in the kinetic equation Eq.([7]), but the microscopic electromagnetic fields are 
fluctuating in a wide dynamical range due to collective plasma efi^ects. The 
coarse grained distribution function of the CR particles / =< / > obeys the 
equation that can be obtained by averaging the microscopic equation Eq.([7]) 
over an ensemble of appropriate short-scale fluctuations 

Here / = / + /', B = B -h b', E = E + E', B =< B >, E =< E > - 

are the averaged flelds, and therefore (b\=0, (e\=0. The ensemble of 

fluctuations can be of external origin or produced by the same population of 
charged particles we only assumed at this point that the collision operator 



/[/,/] = -e E .^ +^ b -Of), (9) 




is a functional of the averaged distribution function / and can be expressed 
through the statistical momenta of the fluctuating field. The collision oper- 
ator describes the momentum and energy exchange between CRs and the 
background plasma and therefore it must be accounted for in the averaged 
governing equations for both the CRs and background plasma. 
The momentum exchange rate is the first moment of Eq. ^ 

p/[/]rf3p = -e ^uIe) + i (^C X b') , (10) 

where n^j., j^^ - are the fluctuating parts of the CR number density and the 
CRs electric current defined by 

''f'd'p, (11) 

{p)f'd'p, (12) 



where v(p) - is the CR particle velocity, and (/ ) = 0. 
Then, by averaging the last term in Eq.®, one can get 



^((5 



ericrU ) X B ) = -(j'^'' - ericru) x B - e (n[.^E,'\ + - /j|,,. x b' 

(13) 



where ricr, jcr - are the averaged CR number density and their electric current, 
jcr = jcr + Jcrj '^cr = «cr + '^cr' Note that Eg.lfTU]) and the last two terms on 
the right hand side of Eg. p^ are coincident. Therefore, we conclude that the 
CR scattering due to the stochastic electromagnetic fields accounted for in 
the kinetic equation Eq. ([8]) by the collision operator must be simultaneously 
included into the equation of motion of the background plasma using Ea. (|13p . 
The averaged induction equation Eq.([5]) can be expressed as 

— = V X (u X B), (14) 

and the averaged equation of motion Eq.([6]) for the background plasma 



du 



(15) 
where pg -is the averaged pressure of background plasma. Note that Eq. p^ 
and Ea.([T5|) is also valid for CRs consisting of electrons and positrons, with 
ricT being the difference between the positron and the electron number den- 
sities, while y^ - the total electric current of the particles. 

In a few cases, namely, for weakly fluctuating magnetic fields or, for strong 
magnetic fiuctuations but at scales smaller than the CR gyroradii, some clo- 
sure procedure s exist to reduce the coUision operator /[/, / ] to /[/] (see e.g. 
lToDtvginlll983l [Bvkov et al.ll2012f ). It is instructive, nevertheless, to derive 
the force density J pl[f](fip for the most simple case of /[/]. The simplest 
form of the collision operator is the relaxation time approximation in the rest 
frame of the background plasma 

m^-Hf-fiso), (16) 

where /iso - is the isotropic part of the momentum distributio n /, and v is 
the co llision frequency due to CR particle- wave interactions fe.g. lBvkov et al.l 
l201ld ). This approach usually implies that the scatterers have no mean (or 
drift) velocity relative to the rest frame of the background plasma. This is not 
always true, if the plasma instabilities that are producing the magnetic field 
fiuctuations are highly anisotropic. However, it can be used to illustrate the 
importance of the momentum exchange between CRs and the background 
plasma. 

Using the parameterisation ly = ail, where fl = — - — , Bq is the mean 

magnetic field, and a - is the CR coUisionality parameter, from Ea. (jl6[) . one 
can obtain 

p/[/]d3p=-^j„.. (17) 

This is the force density in Ea.([T5|). 



3 Instabilities driven by anisotropic CR distributions: 
the kinetic approach 



Consider incompressible modes propagating along the mean homogeneous 
magnetic field Bq in the rest frame of the background plasma. The lin- 
ear dispersion relation can be obtained by the standard perturbation anal- 
ysis of Ea. lfM)) . Ea. ((T5|) and Eq.®, assuming the small perturbations of 
magnetic field b, plasma bulk velocity u and the CR distribution / to be 
ex exp (ikx — itut). The unperturbed anisotropic CR distribution, that is the 
source of the instability free energy, can be represented as 






rN{p) 

47r 



X 



l + 'iP^i + ^ (3a*2 - 1) 



(18) 



where /i = cos 9 , 6 - is the CR particle pitch-angle, ricr- CR number density. 
The multipole moments of the CR angular distribution are parameterized by 
/? (the dipole) and x (the quadrupole). We assume below /3 < 1 and x ^ 1- 
The unperturbed state can be a steady state of a system with CRs where 
both the anisotropy and the spectral distribution N{p) are determined by the 
energy source and sink as well as the magnetic field geometry through the ki- 
netic equation Eq.® with some appropriate boundary conditions. The most 
interesting applicatio n of the formalism is related to diffusive shock acceler- 
ation model (see e.g. Blandford and Eichlerlll987l iMalkov and Drurvll200ll 
iBvkov et al.ll2012 . ISchure et al.ll2012n . In that case the normalized power-law 
CR spectrum is appropriate: 



N{p) = y 



{a - 3) p[,"" 



-3) 



1 



(S)' 



PO<P< Pm, 



(19) 



p" 



where a - is the spectral index, po and Pm - are the minimal and maximal CR 
momenta, respectively. In the DSA applications it is convenient to express 

Us 

the dipole anisotropy parameters through the shock velocity Ug as /3 = — . 

c 
Then dispersion equation has the form: 

I / 47rP77 \' 

iJ- = vl Ik^^k (1 ± ia) ( fcoAo (xq, a;™) H —^^^^Ai (xq, x^) ) - ko 



where Va - 
kcpn 






C I L f r Lli a m Jb 



(20) 
kcp kcpo 



eBo ' 



^0,1 ixo,Xm) = / 0-0,1 {p)N {p)p dp 

Jpo 
oi (P) = 7 / , , ,^ d.n, 



4 J_i 1 =F x/i ± ia 



(21) 
(22) 
(23) 



10 



where the ± signs correspond to the two possible circular polarizations de- 
fined by b = b{ey ± ie^), with the x-axis along the mean field Bq. The func- 
tions Aq^i {xQ,Xm.) are expressed in elementary functions in Appendix A. In 
the collisionless limit a — > the contribution of the pole to the i maginary 
part of Eg . (1221) describes the well known resonant instability (e.g. IZweibell 
I1979L lAmato 2011) , while the real part (the priiicipal p art of the integral) is 
respon sible for the instability discovered by iBelll ( 2004) (see also I Achterbergj 
(|1983D 1. 



The kinetic approach we used here to derive the dispersion equation al- 
lows us to unify the instabilities due to both the dipolc and quadrupole- 
type CR anisotropy. The finite mean free path of the CRs is characterized 
by the coUisionality parameter a. The approach used above allows one to 
study the instabilities driven by the CR anisotropy for arbitrary relations 
between the mode wavelength, the CR mean free path and the CR gyroradii. 
It is instructive to demonstrate the transition between the collisionless case 
(i.e. a = 0), where the CR mean free part is much larger than the mode 
wavelength, and the opposite case with the coUisionality parameter a — >■ 1 
(Bohm's diffusion limit). In the collisionless limit (i.e. a = 0) the instabi l- 
ities due to dipole ty pe ai iisotropy (y = 0) were discussed by iBelll (|2004[) , 
iPelletier et"aLl (|2006l) . and lAmato and Blasil (|2009l) . The firehose instability 



of a highly relativistic plasma without a dipol e anis otropy was discussed 
bv iNoerdlinger and Yuil (|l968h . ISchure and Belli (|201lh derived a dispersion 



equation for the mono-energetic particle distribution instead of the power-law 
distribution in Ea. (|19p used here, and the dipole-type initial anisotropy (i.e. 
X =0). The fireh ose instability of the anisotropic CR pressure with nonzero 
X was studied by iBvkov et al.l ( 2011b[) . 



4 Grov^rth rates of incompressible modes propagating along the 
mean magnetic field 

In Figure [T] we illustrate the growth rates derived from Eq. (PU)) for a par- 
ticular choice of parameters of the CR distribution functions typical for the 
upstream distribution of CRs accelerated by the diffusive acceleration at a 

shock of velocity — — 0.01, with a = 4, and — ^ = 100. The DSA spectrum 

may span many decades, but we choose the two-decade range of the parti- 
cle spectrum to model the instability far upstream of the shock where the 
longwavelength fluctuation amplification is the most efficient. The CR dis- 
tribution function and the CR current normalizations are fixed here by the 

dimensionless parameter fco^'go = 100, where r^o = .To estimate the nor- 

eoQ 

malization of the CR distribution we assumed that about 10% of the shock 

ram pressure is converted into the CR energy. For the CR spectrum of the 

index a = 4 the fraction of CRs above the momentum po is oc Pm/po, while 

Tgo oc po- Therefore the spatial dependence of the key governing parameter 

of the Bell instability fegrgo depends basically on the energy dependent CR 

anisotropy. The bulk of the CRs are confined in the accelerator and therefore 



11 



would have anisotropy about Us/c (apart from the particles at the very end 
of the CR spectrum escaping from the system). 



4.1 Nonresonant shortwavelength instability 

It is instructive to consider the short-scale CR-current driven modes produced 
by Bell's instability as an asymptotic case of the general Eg. (1^111) . f or diff erent 
wavenumber s k in t he coUisionless case a = 0, following iBelll ( 2004[ ) and 
iBvkov et al.l ( 2011br) . In Figure [U we illustrate the growth rate dependence 



on the coUisionality parameter. 

In the wavenumbcr range /coTgo > krg^ > 1, corresponding to the insta- 
bility discovered bv lBelll (J2004l ). the growth of the right hand polarized mode 



(panel a in Figure [T]) is inuch faster than the left hand mode (panel b in 
Figure [IJ . This results in fast helicity production. In the coUisionless limit 
the right hand mode has the growth rate 



76 == Va\/kok~k'^. (24) 

Eq. ((24| follows from Eq. ((20| , neglecting the response of the CR current on the 
magnetic fluctuations, i.e., Ao{xo,Xm) -^ and Ai{xQ,Xm) ^ 0. The weak 
CR-current response is the main cause of the Bell-type instability. Indeed, 
the CR current induces the compensatory reverse current in the background 
plasma and if the current is not responding to a magnetic field variation, then 
the magnetic fluctuation is growing due to the Ampere force. The CR current 
only weakly responds to the magnetic field fluctuations with wavenumbers 
korgo > kvgo > 1, and they grow. From Ea.(|24p one may see that 7^ ~ k^ 
for k -^ ko. 



4.2 The resonant instability 

In the coUisionless case for the wavenumber regimes Xm > 1, but xq < 1, 
the resonant contribution dominates the pole in the integrand in Ea.p2p. 
Therefore, the resonant mode growth can be seen in Figure [T] in the regime 
0.01 < krgQ < 1, where both circular polarization modes are growing with 
the very close rates ex k for a = 4 (compare panels (a) and (b) in Figure [IJ. 
Collisions do not change the mode growth drastically for a < 0.1, but in the 
limit of strong collisions with a = I the left hand mode grows slower than 
the right hand polarized mode. This may also result in helicity production. 



4.3 A nonresonant longwavelcngth instability: the firehose mode 

In the longwavelcngth regime where Xm = ^ 1, within the coUision- 

eBo 

less case, the dispersion relation in Eg. ipO)) can be approximated, following 



12 



Im[uj/{varJ)] 



W 



(a) 
10^ 

10° 

10-^ 

10-' 

10"^ 

10-* 

10" 

(b) 

10^ 

lo" 

10-^ 

10-' 



10"' 









..^■^ 


-^\ 








/ 


/> 




\ 




/ 


/ 






\ 

\ 


/ 


/ 








\ 


// 














\y 








































y 

// 


:'-x 


X 






• 

// 


/ 




\ 


\ 

\ 


i 

/ 

// 

// 


'/ 

/ 








\ 


/ 

/ 













IQ-' 



10-' 



10-' 



10" 



10^ 



10^ 



10' 
kr 



sO 



Fig. 1 The growth rates for the two circularly polarized modes. The right hand 
polarized mode (panel a) and the left hand mode (panel 6) are derived from EQ. (|20p . 
We illustrate the growth rate dependence on the coUisionality parameter a. Dotted 
line corresponds to a=0.01, dashed line - a=0.1, and dot-dashed line - a=l. The 
quadrupole anisotropy is x = 6 (its/c)'. Note that in the bottom panel the dashed 
and dotted lines are very close. 
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iBvkovet al.l(|2011b[) . as 

(25) 



co' = vlk'\lT'f 






Bq (i_ bl 

Pm 



As it follows from Eg. ipS)) . in the regime dominated by the dipole CR anisotropy 
(y — >■ 0) o nly the left-polarized mode is growing with the rate ex fc^ (see 
ISchure and Bell 20 n,). For a finite quadrupole-type CR anisotropy x at small 
enough wavenumbers the modes of both circular polarizations arc growing 
again with the very close rates oc k (see in Figure [J). The instability due to 
the quadrupole-type CR anisotropy corresponds to the well known firehose 
instability in a plasma with anisotropic pressure. Indeed, the CR pressure 
anisotropy derived from the CR distribution Ea.([T^ is 

Pf - PZ = ^xP'^ (26) 

where 

1 1"°° 

P^'' = x"cr/ v{p)N{p)p'dp. (27) 

-J Jo 

The dispersion relation for the modes produced by only the quadrupole-type 
anisotropy of CR distribution can be obtained from Eg . ([251) if one neglects 
the dipole-type contribution Xm -> 0. Then, it is reduced to the standard 
hydrodynamic dispersion relation of the firehose instability 




(28) 

where P\\ — P± - i s the pressure anisotropy along the mean magnetic field 
direc tion (see, e.g.. lBlandford and Eichlerlll987LllYeumann and BaumiohannI 
Il997[ ). The dispersion relation Ea. ipS)) is justified for the modes with the 
wavenumbers above the CR ion gyroradii. The dependence of the growth 
rates of the firehouse instability on the collisionality parameter can be seen 
in Figure [1] It should be noted that the growth rates of the firehose modes 
of both polarizations in the regime krgo < 1 are declining functions of the 
collisionality parameter. Their growth rates would be equal in the case of lack 
of the CR current. Contrary, the growth rates of the current driven modes 
are different for the two polarizations. The growth rate of the right hand 
pola rized CR-c urrcnt driven mode is sensitive to the collisionality parameter 
(see ISchure and Bcll,2011.1 . 



4.4 A nonrcsonant long-wavelength instability: the cosmic-ray current 
driven dynamo 

Bell's instability results in the fast growth of short-scale modes with wave- 
lengths shorter than the gyroradius of the cosmic-ray particles and in the 
presence of CR-current it may produce strong short-scale turbulence (e.g. 
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Fig. 2 The growth rates of the longwavelength modes of two circular polarizations. 
The right hand polarized mode (panel a) and the left hand mode (panel 6) are prop- 
agating along the mean magnetic field as function of the wav enumber. The dotted 
line curves are derived from the dispersion equation Ea. H3ip for the collisionality 
parameter a — 0.1, the dimensionless r.m.s. amplitude of Bell's turbulence Nb = 1, 
and the mixing parameter ^ = 3. The dashed curves given for comparison are the 
growth rates derived from Ea. (|20p which are shown in Figure [T] 
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Bell and Lucekl200lllBelJ2004IZirakashvili and Ptuskiij2008llZirakashvili et alJ 



2008L iReville et al.ll2008l IVladimirov et alJl2009L iRogachevskii et all 12011) 



Moreover, the shortscale turbulence is helical, and at the wavenumbers be- 
low 2 ko its kinetic energy density dominates over the magnetic e nergy density 
making a favorable condition for a pure a-dynamo effect (see iBvkov et al.1 
[2011c) . The strong short-scale turbulence influe nces the background plasma 
dynamics on scales larger than the CR gvroradii. lBvkov et alJ ()2011ct ) derived 



the mean field dynamic equations averaged over the ensemble of short-scale 
motions for plasma systems with CR-current. The averaged equation of mo- 
tion can be presented as 

i9V 1 1 

— + (VV)V = --VPg - ((uV)u) + _ ((V X b) X b) + 

+ -^((V X B) X B) - - — ((j" - en„.V) x B) - f pl[f]d^p, (29) 
inp cp J 

where V is the mean velocity of the plasma. The magnetic induction equation 
for the mean magnetic field B reads 

— = cVxf + Vx(VxB) + i^,„AB. (30) 

Here £ ~ (u x b) is the average turbulent electromotive force and i/m is 
the magnetic diffusivity. The averaged equations Ea.ip^ and Eg. ([5(1)) are 
designed to be applied to the dynamics of modes with scales larger than r^o , 
i.e., CR particles arc magnetized on these scales. 

The ponderomotive forces ((uV)u) and -^ ((V x b) x b) in Ea. ([29[) de- 
scribe the momentum exchange of the background plasma with the Bell 
mode turbulence. The averaged turbulent electromotive force £ results in 
the magnetic induction evolution. It is important that in the case under con- 
sideration the ponderomotive forces in Ea. ([^^ depend on the CR current 
through the Bell mode turbulence moments. To e xpress the electrom otive 
and ponderomotive forces through the CR current [Bvkov et al.l ( 2011c[ ) fol- 



lowed the mean field closu re procedure similar to the approach proposed by 
Blackman and Fieldl ([2002[ ) in the dynamo theory (see for a review [Brandenburg| 



2009af ). The closure procedure is introduced by the parameter Tcor- The cor- 
relation time Tcor which is the relaxation time of triple correlations and is 
approximately equal to the turnover time of the Bell turbulence. The de- 
pendence of the electromotive force and the ponderomotive force on the CR 
current is determined by the kinetic coefficients at and Kt, correspondingly. 
The kinetic coefficients are determined by the r.m.s. amplitude of Bell's tur- 
bulence (6^) and Tcor. The short scale turbulence produced by the Bell mode 
instability is helical and therefore there is also a contribution to the electro- 
motive force oc atB resulting in the a-dynamo effect. Then, the dispersion 
equation for the modes of wavelengths longer than r^o in a plasma with 
anisotropic relativistic CRs can be derived from Eg. ([^5]) and Eg. ([50)) by the 
standard linear perturbation analysis: 



2 I 2 2 -^ ., "t 

uj — k v^^ Lotk 

Airp 



1 /, . , , AnericrX , , \\ 3, 

- I fcoAo (a;o, a;,„) H Ai [xq, Xm) I + -kq 
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i>0 



+ iakvl ( koAo (xq, x™) H —^^—Ai (xq, x™) ) = 0. (31) 



The dispersion relation Eq.ljSTj) was derived for the systems where the un- 
perturbed CR-current is directed along the unperturbed magnetic field, and 
the short scale turbulence consists of Bell's modes . It is convenient to intro- 

duce two dimensionless parameters Nb ~ — — — Bell's turbulence r.m.s. 

-So 
amplitude, and the dimensionless mixing length ^, instead of the correla- 
tion time Tcoi-. The mixing length is defined here as 27r^/fco = Tcor\/ {v'^) ~ 

Tcov^/UJWIWp)- Then at w (6|)Tcor ~ Sir^y^NBVakQ^p and Kt = ttNeBq. 

In Figure [5] we illustrate the long wavelength mode growth derived from 
Ea. (|5T|) for ^ = 3. The corresponding mixing length is close to the scale of the 
maximal growth rate of the short scale Bell's instability. The a-dynamo ef- 
fect dominates the growth rate of a polarized mode shown in Figure [2] (panel 
b) in the intermediate wavenumber regime a < kvgo < 1. One should have 
in mind that in the case of Bohm's CR diffusion a ~ 1 and therefore the 
intermediate wavenumber regime is rather limited. It should be noted that 
the helicity of the unstable, long- wavelength mode studied above is opposite 
to that of the short-scale Bell mode. This provides, at least in principle, the 
possibility of balancing the global helicity of the system by combining short 
and long-wavelength modes. Care must be taken however, since numerical 
models indicate a high saturation amplitude of the Bell mode making a non- 
linear analysis necessary to address the helicity balance issue. We will discuss 
some nonlinear simulations below in fj5l 



4.5 The cosmic-ray current driven instability in the hydrodynamic regime 

The nonresonant modes in a hydrodynamic regime, where the wavelength is 
longer than t he mean free path, i .e., krgo < a, are unstable, as it follows from 
Ea. (|3T|) (see iBvkov et al.|l201ld . for details). Both circular polarizations in 
panels a and b in Figure [5] grow with the same rate given by 



7 ~ \l —^Vkkoava- (32) 

The transition from the intermediate wavenumber regime a < fcr^o < 1, 
dominated by the dynamo effect discussed in ij4.41 where the mode growth 
rate can be approximated by 



7 



ATT^NBVak, (33) 



to the hydrodynamical regime with fcr^o < a where the growth rate is oc fc^' ^ 
according to Eq. ([5^ , is clearly seen in Figure [5] (panel b) . Note that for the 
mode polarization shown with the dotted line in the panel a of Figure [2l no 
dynamo-type instability occurs, but the hydrodynamical regime instability 
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is present. This mode grow fast in the short wavelength regime krgQ > 1 due 
to BeU's instability. 

The effect of the short-scale turbulence on the hydrodynamic regime insta- 
bility enters Eq. (PT|) through the turbulent coefficient nt/Bo. The turbulent 
ponderomotive force is large enough in both the intermediate and hydrody- 
naniical regimes, and the CR current response in the long-wavelength regime 
can no longer be neglected. The current cannot be treated as a fixed exter- 
nal parameter, as is normally done for the short-scale Bell instability, and 
therefore the MHD models of the Bell tur b ulence that assume a constant 
CR-current (see e.g. Bell and Lucekl 2001 . Zirakashvili an d Ptus kin 20081 



Zirakashvili et al.l2008llReville et al.l2008llVladimirov et al..2009K,Rogachevskii et al.l 



2012f ) cannot be directly applied to the nonlinear models of the longwave 



length instabilities discussed above. Particle- in-cell simulations with very lim- 
ited dynamical range performed by iRiouelme and Spitkovskvl (J2009l l201Gl ) 
indicate the importance of the CR backreaction effect on the CR-driven in- 
stabilities. Therefore the nonlinear dynamics of the long- wave CR-driven tur- 
bulence in a wide dynamical range remains to be investigated. In the next 
section we illustrate the nonlinear evolution of the short scale turbulence 
driven by a fixed CR current, using high resolution MHD simulations. 



5 Numerical solutions of the Bell dynamo instability 

Significant insights have been possible through high-resolution direct numer- 
ical simulations (DNS) and large eddy simulations (LES) of the Bell insta- 
bility and its subsequent saturation. In this section we describe some of the 
main results and, in particular, the connection with the dynamo instability. 
The simulations have been carried out in a Cartesian domain of size L"^, 
so the smallest wavenumber in that domain is fci ~ 27t/L. The system is 
characterized by the non-dimensional parameter 

47r j" , , 

C Ki±Jo 

In the ideal case (yu = 0), the Bell instability is excited when J > \ and 
the normalized wavenumber of the fastest growing mode is fc/fci = J jl. 

The normalized growth rate of this fastest growing mode is jh/^Aoki = 

j7/2. In Figure [3] we reproduce the results of numerical simulati o ns of iBelll 
(|2004[ ) for J^ = 2 using 128^ mesh points and IZirakashviU et al.l (|2008[) for 



j/ = 16 using 256"^ mesh points. These simulations confirmed the analytically 
expected linear growth rates. Interestingly, the saturation of the instability 
was never perfect. Instead, the ma gnetic field still continued to grow at a 
slow rate. iRoeachevskii et al.l ( 2012t) have argued that this slow growth after 



the end of the exponential growth phase of the instability is the result of 
a mean-field a effect. The purpose of this section is to elaborate on this 

possibility. 

We begin by discussing first the recent DNS of lRogachevskii et al.l (|2012f ) 
for J' = 80 and J = 800 at a resolution of 512'^ mesh points and discuss 
also new results for J' = 800 at a resolution of 1024"^ mesh points. In all 
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cases, explicit viscosity v and magnetic diffusivity v^ are used, so the fastest 
growing modes in those cases have somewhat smaller wavenumbers than in 
the ideal case. This is quantified by the Lundquist number Lu = Valvyiki and 
the id eal case corresponds then to Lu ^- oo. For example. iRogachevskii et alj 
( 2012[ ) used Lu = 80, in which case the fastest growing mode has fc^/fci sa 21 
for J = m while for J = 800 it has /cjfci « 63. The DNS show that most of 
the power is at somewhat larger wavenumbers; see Figured where we show 
magnetic energy spectra for both cases. 

In Figure [5] we show the temporal evolution of spectral magnetic en- 
ergy Em and the spectral kinetic energy Ek at selected wavenumbers. These 
curves show an exponential growth at early times, followed by a slower growth 
at later times. At the wavenumbers of the Bell mode, the growth rate from 
linear theory is reproduced. At smaller wavenumbers, the growth is at first 
slower, and then it is even faster than the growth rate of the Bell mode. This is 
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Fig. 3 Numerical solutions of the Bell instability for J = 2 using 128^ mesh points 
(IBelll[200l . left hand side) and J = V6 using 256^ mesh points (jZirakashvili et aP 
120081 . right hand side) . Note the continued growth of the magnetic field at the end 
of the linear growth phase at t « 10 on the left and i ~ 1 in the right. Courtesy of 
Tony Bell (left panel) and Vladimir Zirakashvili (right panel). 
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Fig. 4 Time evolution of Euik.t) for J" = 80 (left) and J = 800 (right) at 
resolutions 512'' and 1024'^, respectively. The solid lines refer to the initial spectra 
proportional to fc'' for small values of k and the red and blue lines represent the 
last instant of Eu and Ek^ respectively. The straight lines show the fc* and k~^'^ 
power laws. 
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Fig. 5 Time evolution of Em ki/v'^g for ^ = 80 (left) at wavenumbers fc/fci = 1 
(solid line), 5 (dotted), and 21 (dashed) and J = 800 (right) at wavenumbers 
k/ki — 1 (solid line), 10 (dotted), and 63 (dashed). The short straight lines show 
the growth of the energies for the Bell (dashed) and dynamo (solid) instabilities. 



a consequence of mode coupling ( Rogachevskii at ani2012f ). Comparing with 
Figure |4l we see that after some time a k^ energy spectrum is established. 
Such an energy spectrum is also known as Batchelor spec trum and can be de- 
rived under the constraints of solenoidality and causality (JDurrer and Caprinil 
|2003[ ). When the fc^ spectrum is established, the growth of spectral energy at 
small wavenumbers is no longer described by linear theory, but follows the 
growth of the Bell mode. 

In Figure |6] we show visualizations of B^/Bq on the periphery of the 
computational domain for J' — 80 using 512'^ mesh points and J^ — 800 
using 1024^ mesh points at two different times. One clearly sees that at early 
times, the magnetic field shows a layered structure with a high wavenumber 
in the z direction. At later times, the magnetic field breaks up and becomes 
turbulent. In both cases, larger scale structures develop, as one also sees from 
the energy spectra in Figure SI 

It should be pointed out that, owing to the persistent growth of magnetic 
and kinetic energy, the Reynolds numbers grow eventually beyond the limit 
of what can be resolved at a given resolution. Unlike some of the earlier 
LES, where numerical ef fective viscosity and diffus ivity keep the small scales 
resolved, in the DNS of iRoeachevskii et al.l (J2012t) this is not the case and 
the numerical code (in this case the Pencil CodeQ) eventually 'crashes'. 
The point when this happens can be delayed by using higher resolution. This 
is why we show here the results for J' = 800 at a resolution of 1024^ mesh 
points, where the simulation can be carried out for about 0.126 Alfven times, 
co mpared to only 0.09 Alfvc n times at a resolution of 512"^ mesh points used 
in lRoeachevskii et al.l ( 20121) . Remeshing the 1024'^ run to 2048'^ mesh points, 
we were able to continue until 0.142 Alfven times, after which we were unable 
to continue the run due to a disk problem. 

The Bell instability is driven by the simultaneous presence of an exter- 
nal magnetic field Bq and an external current jcr, giving therefore rise to a 
pseudo-scalar jcr • Bq; here, Bq is an axial vector while jcr is a polar vector. 
In stellar magnetism, the presence of a pseudo-scalar is caused by rotation 
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Fig. 6 Visualization of Bx/Bo on the periphery of the computational domain 
^ = 80 using 512^ mesh points (upper row) and J = 800 using 1024^ mesh points 
(lower row) with Lundquist number Lu = 80 in both cases. 



n (an axial vector) and gravity g (a polar vector) . This property is generally 
held responsible for the production of magnetic fields by what is known as 
the a effect. As explained in Sect. 4.4, the a effect denotes the presence of 
a tensorial connection between a mean electromotive force £ = u x b and a 
mean magnetic field via 



aijBj + rjijkBj^k 



(35) 



where higher order derivatives (indicated by commas) of the mean magnetic 
field are also present. If the tensors a^- and rjijk were isotropic and the evo- 
lution characterized by just two quantities, a = bijOiijji and ?/t = ^ijkflijk/^^ 
the growth of the mean magnetic field would occur at the rate 



7dyn 



ak — rjTk , 



(36) 



where rjT = i^t + vm is the total (turbulent plus microphysical) magnetic 
diffusivity and the fastest growth occurs at wavenumber k = a/2r]T with the 



growth rate 7n 



a^/Arn 
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In stellar dynamos, where the magnetic Reynolds number is very large, 
the actual growth is dominated by small-scale dynamo action, so Ea. (|36p 
is in practice not obeyed, unless the small-s cale dynamo is not excited, for 
example at low magnetic Prandtl numbers ( Brandenburg! [2009bl ) . However, 



in the present case the magnetic energy spectra show that at late times, 
magnetic power moves gradually to larger scales. This is why we now ask 
wh ether this can be explained by the a effect. 

iRogachevskii et aD ( 20121 ) have shown that in the case of jcr and Bq 



pointing in the z direction, the large-scale mean magnetic field is a function 
of X and y and can be written in terms of two scalar functions ^ii (x, y, t) and 

-B|| (x, y, t) with 

B(x,2/,i) = Vx(zA||)+£Bj|, (37) 

where z = (0, 0, 1) is the unit vector in the z direction. These functions obey 
the mean field equations 

5:4,1 /dt = a^B II -f- 77^ V^ A|| , (38) 

dB\\ /dt = as J|| + rys V^:B|| , (39) 

where Jii — — V^^ii is the xy dependent part of the mean current density in 
the z direction. We consider a homogeneous system, so the coefficients a^i, 
as, rjA, and tjb are constant and we can seek solutions of a form proportional 
to exp(Af -I- ik • x). In this case, the dynamo growth rate is still described by 
Ea. ([55)) . provided we substitute 



a 



=« = (aAaB -1-62^2)1/2 ^j^j TJt-^Tjf^{TjA+ljB)/2, (40) 



where e,, = {tja ^?7b)/2 quantifies the anisotropy of the turbulent diffusivity. 
To determin e these coefficients from the DNS, we use the so-called test- 
field method of ISchrinner et al.l ( 2005^ . which was originally used in spher- 



ic al coordinates. The i mplementation in Cartesian coordinat e s is d escribed 
in iBrandenburd ( 2005h and especially in [Brandenburg et ahl (|2Q12|) , where 



the mean magnetic field was allowed to depen d on all three s patia, l coordi- 
nates, and not just on one, as was assumed in IBrandenburd (l2005f ). Under 
the assumption that the turbulence is governed by only one preferred di- 
rection, which is here the case, the number of coefficients reduces to 9, and 
homogeneity reduces this number further to 5, so in the present case we have 

£ = a±B_L + a||B|| - /3_lJ± + /3|| J|| - ^z x K_l, (41) 

where J = V x B characterizes the antisymmetric part of the magnetic deriva- 
tive tensor and Ki = (Bji + Bi,i)zi/2 the sym metric part. We have followed 
here the notation of iBrandenburg et al.l (|2012l ). except that there the two a 



coefficients were defined with the opposite sign. Comparing with the coeffi- 
cients used in Eas. ([55)) and ([M|. we find that a^ = a||, as = ci±, Va = l3\\, 
and r]B = /3_l — m/2- In Figure [7] we show the time dependence of the various 
parameter combinations. In the early kinematic phase (tvAoki < 0.08), the 
root mean square velocity, Umis, as well as a\\ and a± grow exponentially. At 
later times, ay continues to grow, while a± remains small and approximately 
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Fig. 7 Time evolution of the model parameters for J = 800 and Lu = 80 using 
1024^ mesh points, (a) Exponential growth and subsequent near-saturation of Urms, 
«ll, and a± (all normalized by wao) in linear-logarithmic representation, (b) Evo- 
lution of a|| and a± (normalized by «ao) in double linear representation, showing 
that a± is much smaller than Q||. (c) Evolution of /3||, /3_l, and fj, (normalized by 
VAo/ki). (d) Evolution of ?7t and e^ (normalized by VAo/ki). (e) Evolution of Ca 
(negative values are shown as dotted lines), and (f) growth of the fasted growing 
mode. 



constant. The other turbulent transport coefficients also grow exponentially 
in the kinematic phase, and at later times /3||, rjt, and e,, continue to grow, 
while I3± and fi remain small and can even become negative. The resulting 
effective dynamo number, which is proportional to the product anax, reaches 
values well above the critical value of unity. The estimated and actual growth 
rates agree roughly and have a value of around 10 in units of ijtki. 
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6 Instabilities driven by the nearly isotropic CR distributions 

In many astrophysical objects the CR mean free path due to the particle 
scattering by magnetic fluctuations carried by the background plasmas is 
below the characteristic scale sizes of the plasma flow. In that case the angu- 
lar distribution of the CRs is nearly isotropic with a small anisotropic part 
(i.e. both /? <C 1 and x *S 1 in Ea. ([T5)) ). Then one can use the diffusion 
approximation that assumes 



^^^(^'P)-i 



7V^'-(r,p) + — pJ-(r,p) 
vp 



where the diffusive current of CRs is 



KaflVsA^'^ 



pdN" 



3 dp 



(42) 



(43) 



Kap is the momentum-dependent CR diffusion tensor. Then the kinetic equa- 
tion Eq.® reduces to the advection-diffusion equation for the isotropic part 
of CR distribution Ar'=''(r,p,i) 



dN" 
~dt 



OL^a.13 



VpN""" - (uV) N" 



pdN" 



-Vu, 



(44) 



wher e u(r, t) is the bulk velocity of the background plasma (see, e.g.. lToDtvginl 
Il983f ). It is assumed here for simplicity that the scattcrers are carried with 
the plasma bulk veloci ty, though it is possible to account for the scatterers 
drift velocity (see, e.g.. lSkilling|ll975l) . The advantage of this approach is that 
it is valid for collision operators /[/] more general than just the simple relax- 
ation time approximation given by Ea. (jl6p . In the diffusion approximation 
the exact form of the collision operator determines the form of the diffusion 
tensor and its momentum dependence. Therefore, the results obtained within 
the diffusion approximation are valid for different collision operators. 

To explore the effect of CRs on the background plasma one should cal- 
culate the first moment of the kinetic equation Eq.® for CRs that is the 
momentum exchange rate between the CRs and the background plasma: 



dVc. 



V„P- + V^TIL 



^{r-en,,u)xB + JpI[f]d^p 



where P'^'' is the CR pressure, the CR momentum density 

P(r, t) = jpfd'p, 
and the reduced CR momentum flux density W^n is defined by 



, (45) 



(46) 



.;.^/../.'.-..., 



(47) 
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In the diffusion approximation for the steady state (e.g., in the shock rest 
frame) the first and the third terms in the left hand side of Eq. (|^5|) are small 
and then Eq. p3| can be reduced to 



| + (uV)u)^-V(,, + P-) + ^. 



p{— + (uV)u ) = -V {pg + P-) + — (V X B) X B. (48) 



The equation can be applied to longwavelength perturbations. It should be 
supplied with the continuity equation: 

|^ + V(pu) = 0, (49) 

the energy equations for the background plasma: 

^ + (uV)pg + 79PsVu = 0, (50) 

the MHD induction equation 

— = V X (u X B) , VB = 0, (51) 

and the equation for CR-pressure variations 

gpcr 

(uV)P^'■+7,,P^'^Vu = V„K„^V^P^^ (52) 



dt 



where Kap is the CR diffusion tensor averaged over the CR distribution func- 
tion, 7g and 7c,- - are the adiabatic indexes of the plasma and CRs, respec- 
tively. 



7 Acoustic instability driven by the CR pressure gradient 



It was found by 'DrurvV 1984!) jDorfi and Drurvl(|l985[ ). lDrurv and Fallel(| 19861 ) 



iDrurv and Downes (2012 ) that the force density in Ea.(|^5| associated with 
the CR pressure gradient that does not depend on the density of the back- 
ground plasma results in a speci fic instabi li ty. Th e effe ct of ma g netic fi eld on 
the instability was studied by iBerezhkol (|1986[) and IChalovl ([19813) • The 



analytical study of the instability can be performed for the modes with 
the wavenumbers below the scale size of the CR pressure gradient L ~ 
per ^|y per I ^^ |-j^g gene ric case of the d i ffusiv e shock a c celerat ion L ~ 
(c/us) X rg/a. Following iDrurv and Falli (|1986[) . IChalovl (|l988bl ) for the 



wavenumber range kL > 1, but krg/a < 1 the mode growth and damping 
can be derived from the continuity equation for the wave action. 

The mode growth rate P in the simplified geometry where the CR pressure 
gradient is directed along the unperturbed magnetic field was derived using 
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a standard linear analysis of Eqs. 
the following expression 



r 



in -[52]) bv lChalovl (|l988a[) . who obtained 



ry per 






2vl 



{v1 + vl) 



Po Ko||A:| + Ko_Lfci vl 



± 



± 



PoVm k 



"JKOll/ 



"aW 



KQ\\^ + Ko±k\ vl 



(53) 



Here Vg is the sound speed of the background plasma, P^^ is the unperturbed 
CR pressure, VPq'' is the gradient of the unperturbed CR pressure, fc|| and 
k± are the components of the mode wavevector parallel and transverse to the 
unperturbed magnetic field, respectively, and Ko||, ko-L are the components of 
the averaged CR diffusion tensor. It is assumed that the CR diffusion tensor 



components scale with the background plasma density as n\\.i_ 
phase velocity of the mode is 



The 



vt± 



(w2 + vlY 



kl 



(54) 



The first term in Ea.(|53p is the wave da mping rate due to the irreversible 
stochastic Fermi II C R acceleration effect ( Achterbergll979llBvkov and Toptvghinl 
1X9791 . iPtuskinI Il98lh , while the second and the third terms describe the 
growth/damping of th e modes due to the acoustic instability studied by 
iDrurv and Falld ( 1986() . A more general treatment with a n arbitra r y direc - 
tion of the unperturbed magnetic field was performed bv IChalovl ( 1988b() . 
He accounted for the response of the CR diffusion tensor to both the density 
and magnetic field variations and found that the latter does not change the 
character of the angular dependence of the growth rate significantly. A sim- 
ilar angular dependence of the long-wave mode growth rate d ue to the CR 
current driven instability (discussed above in M.Sp was found bv lBvkov et al.l 
(|2nild) . 

In the space plasma with the modest level of the magnetic field fluctu- 
ations the local CR diffusion is anisotropic. For magnetized CR particles 
(a ^ 1) the diffusion parallel to the mean magnetic field dominates over the 
CR diffusion transverse to the mean field, i.e., k,q\\ 3> kqx- The growth rate 
of the acoustic instability in the anisotropic system is maximal for the modes 

propagating nearly transverse to the mean magnetic field {d — > — ). Here d 

is the angle between the mode wavevector and the the mean magnetic field. 
The angular dependence of the growth rate Eq. ([5^ can be approximated 

by 

Go(^) - ,...,'":',._.. , (55) 



where we used 



COS??, 



cos2 t? -I- Sii sin^ d ' 

sin^. The anisotro py of th e CR diffusion is 



determined by the CR particle magnetization (e.g.. lToptvginiil983i l. that is 
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Fig. 8 The ch ara cteristic angular dependence of the growth rate of the acoustic 
instability Ea. H53|l for a =0.3 (the dashed curve), a=0.2 (the dotted curve), a=0.1 
(the dot-dashed curve), o=0.05 (the rare dot curve). 



the inverse collisionality parameter a, and therefore, ex a . The maximal 

growth rate is therefore achieved for the mode propagating at cosiJmax = a, 
where Gmaxf'^^max) — -Z-- The angular dependence of the growth rate of the 

acoustic instability is illustrated in Figure[5]for various values of collisionality 
parameter a. 



The linear perturbation analysis discussed above is based on the diffusion 
approximation of the CR dynamics in Eqs. (|^ - [5^ and, therefore, it is valid 
for the modes of the wavenumbers above the mean free path of the CRs. 
A numerical model of the acoustic instabi l ity in the nonlinear regime was 
performed recently by iDrurv and Downed ( 2012I) . who found a significant 
amplification of magnetic field. The authors assumed a fixed CR diffusion 
gradient with no response of the CR pressure to the fluctuations, that may 
affect the model results. 
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8 Self-confinement of CRs near their acceleration sites 

Apart from being a central issue for the aceeleration in SNR shocks, the CR- 
driven instabihties are fast becoming an integral part of CR escape models. 
One common difficulty with the observational verification of the proton es- 
cape is that, in contrast to electrons, they likely remain invisible until they 
reach some dense material in SNR surroundings. Only there generate they 
enough 7r° mesons in collisions with other protons and the mesons in turn 
decay into gamma photons which may be detected. Not surprisingly, the es- 
cape of CRs from an SNR is a hot topic of today research in gamma-ray 
astronomy. 

The backbone of the DSA is a self-confinement of accelerated particles by 
scattering off various magnetic perturbations that particles drive by them- 
selves while streaming ahead of the shock. Most important of them were 
discussed at some length in this review. Logically, this process should also 
control the ensuing propagation of CRs, before their density drops below the 
wave instability threshold. Strictly speaking the CR release (escape) from 
the accelerator should be treated together with the acceleration, as it does 
not occur at once for all the particles. But this would be a combination of 
two difficult enough problems and most of the progress in CR escape was 
made by considering it separately from acceleration. 

Remarkably, even within this limited approach, and under rather loose 
formulation of the problem, no consensus on the escape mechanism has been 
reached so far; the dividing lines seem to run across the following issues: 
(i) does the escape occur isotropically or along the local magnetic field? 
(ii) does the scattering by the background MHD turbulence control the CR 
propagation alone or self-excited waves need to be included? (iii) if so, is a 
quasilinear saturation of self-excited waves sufficient or nonlinear processes 
of wave damping are crucial to the particle propagation? (iv) if they are, 
which particular mechanism(s) should be employed? 

Starting with (i-ii) we note that most of the early models, and some of the 
recent ones that target specific remnants, assume isotropic CR propagation 
from a point source impeded onl y by the backgr ound turbulence (one may 
call them test particle models, e .g. Aharon ian an d Atovanlll996llGabici et al.l 



l2009l lEllison and Bvkovl 120111 lGabicill201ll). It should b e note d, however 



that e.g., iRosner and Bodol (1 1996) an d Nava and Gabicil (|2013[) adopted a 



field aligned propagation while iDrurvl (l201l[ ) included the finite radius of a 



SNR shock in the CR escape description. Given the topic of the present short 
review, however, we focus in this section on models that explicitly include 
the self-excited waves. Brief revie ws of o t her a spect s of CR propag ation in 
the galaxy were given recently by iGabicil ( 20 111) and iPtuskinI ( 2012 ) . 



The role of self-confinement effects in the CR escape, their subsequent 
propagation and how these phenomena are treated in different models, can 
be best demonstrated by writing the following equations that self-consistcntly 
describe the CR diffusion and wave generation 

d „ . ■. d KB dPcR , . 

^Pcr(p) = ^-^^ (56) 
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d J dPcR . . 

Here Va is the Alfven velocity, kb is the CR diffusion coefficient in Bohm 
regime, kb = crg/Z, and the time derivative is taken along the characteristics 
of unstable Alfven waves, forward propagating along the field (z- direction): 

d d d 

Ea. (|56p above is essentially a well-known convection-diffusion equation, writ- 
ten for the dimensionless CR partial pressure PcR instead of their distribution 
function / {p,t). We normalized it to the magnetic energy density pv\l2: 

PcR = ^AV/, (59) 

3 pvi 

where v and p are the CR speed and momentum, and p- the plasma density. 
The total CR pressure is normalized to dlnp, similarly to the wave energy 
density /: 



Stt 



^Jl{k)d\nk = ^Jlip)d\np 



Ea. (j57p is a wave kinetic equation in which the energy transferred to the 
waves equals the total work don e by the par ticles, {u + v^) VPcR, less the 
work done on the fluid, uVPcR ( Drurvl[l983[ ) (we neglect the bulk flow ve- 



locity u, here and in Eg . ([55)1 assuming that the active phase of acceleration 
ended by this time) . The above interpretation of the wave generation indi- 
cates that it operates in a maximum efficiency regime. A formal quasilinear 
derivation of this equation assumes that the particle momentum p is related 
to the wave number k by the 'sharpened' resonance condition kp = e Bp/c in- 
stead of the conventional cyclotron resonance condition fcpy — cBq/c (ISkillingj 
Il975[ ). (note that here k = fc||). We assume that dPcR/dz < at all times, so 
that only the forward propagating waves are unstable. The latter inequality 
is ensured by the formulation of initial value problem symmetric with respect 
to z = 0, so we consider the CR escape into the half-space z > with the 
boundary condition dPcn/dz = at z = 0. 

Papers on CR self-confinement discussed below use equations that are 
largely similar to Eas. (|56ll57)) but different assumptions are made regarding 
geometry of particle escape from the source (see (i) above), the character 
and strength of wave dainping F (iv) , and the role of quasilinear wave sat- 
uration (iii) . iFuiita et al.l (|201l[ ) and lYan et ahl (|2012[ ) utilize the isotropic 



escape models (in t his c ase d/dz shou l d be r eplaced by d/dr, etc.) while 
iPtuskin et al.l (J2008[) and iMalkov et al.l (l2012al ) assume that part icles prop- 



agate predominantly along the local large-scale field. Note that lYan et all 
( 2012[ ) considered the escape from an active accelerator (in Ea. ((55)) . one 



should include the fiow bulk velocity, Va ^ Va + u in this case) and, in addi- 
tion, they introduce a step- wise increase in CR diffusivity at a certain particle 
momentum above which particles escape the accelerator. These assumptions 
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Fig. 9 Spatial distribution of CR partial pressure (as a function of ( — z/y/avat, 
multiplied by vj ^fatj kb) shown for integrated values of this quantity YI — 3.6; 
6.7; 10.1 and for for the background diffusivity Dism ~ 10*. Exact analytic solu- 
tions are shown with the solid lines while the interpolations given by Equation (|65p 
are shown with the dashed lines. For comparison, a formal test particle solution 
for i7 = 10.1 is also shown with the dot-dashed line. Note the three characteristic 
zones of the CR confinement: the innermost flat top core, the scale invariant (1/C) 
pedestal, and the exponential decay zone. 



make it difficult to compare their re sults with those of the remaining three 
papers. In these. iFuiita et aD ( 201l|) presented the results of numerical in- 
tegration of Eqs. (|561l57)) (in a spherical symmetry) with neglected damping 
term F. The results indicate a considerable delay of diffusion from the source 
due to a self-confinement. 

However, in the regions where magnetic perturbations are weak, i.e. / <C 
1, the field aligned CR transport is appropriate, as the perpendicular diffusion 
is suppressed, k± ~ ^'^^W ^ ''^H — I'^b/ I ■ Taking into account the condition 
^ISM ^ 1; such regime appears inevitable far away from the source and at 
late times when particles are spread over a large volume and the waves are 
driven only weakly. At earlier times and close to the region of the initial 
localization of CRs, an estimate k,±_ ~ km ~ kb appears to be adequate. 
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Both analytical models by iPtuskin et all (|2008[) and iMalkov et all (|2012a[) . 



however, do not embrace the general case and rely on the assumption k± <C 
Kii thus considering a field-aligned escape. At the same time, they are different 
in further simplifications made, that lead to rather different results, both 
qu antitatively an d qual itatively. 

iPtuskin et all (|2008f ) neglect dl/dt on the l.h.s of Eq.dSTj) thus balancing 



the driving term with the damping term on its r.h.s and assume a Kolmogorov 
dissipation for F, 

r = kVa^I/ {2CKf'^ (60) 

with Ck ~ 3.6 and k ~ l/?'g {p) being the resonant wave number. Therefore, 
only one equation ((55)) needs to be solved which lead to the following self- 
similar solution (in notations and normalization used in Eqs. (|56ll57p ) 

4.3-3/2 
PcH = , (61) 

t'^/^^(j + {kzf /t'^ 

where the dimensionless time t' = [nsk'^/^CK) t, cr = T* (1/4) j-K'^Z^rf, 
and r is the gamma function. The single important parameter this solution 
depends on is the integrated (along the field line) CR partial pressure 



r] = 2k PcRdz (62) 



Therefore, the CR density decays at the source as ex i^^/^ and the flat- 
topped, self-confined part of the CR distribution spreads as z oc i^/^, both 
pointing at the superdiffusive CR transport. The reason is clearly in a very 
strong wave damping due to the Kolmogorov dissipation. For the same rea- 
son this solution does not recover the test particle asymptotic result Pcr oc 
t~-'^/^exp (— z^/4Z?isMi), physically expected in z,t — ?> 00 limit in the inter- 
stellar medium with the background diffusion coeffici ent Dism- 

An alternative choice of damping mechanism is the lColdreich and Sridhan 
(|l997[ ) MHD spectrum, which see ms to be more appropriate in / < 1 regime 



unde r not too strong M HD cascade ( Farmer and Goldreich|[2004ilBeresnvak and LazarianI 



l2008llYan et al.ll2012[ ). The damping rate in this case is 



r = ^-\!j^ (63) 

where L is the outer scale of turbulence which may be as large as lOOpc. Not 
only is this damping orders of magnitude (roughly a factor yJvgjL ) lower 
than the Kolmogorov one but, as it does not depend on / and can be con- 
sidered as coordinate independent, it allows for the following ('quasilinear') 
integral of the system of Eas. ([55| and ([5i 



Pc„(M) = Pc„.(.')-^|.n^ (64) 



Wa dz /o {z') 
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Here Pqro (z) and Iq (z) are the initial distributions of the CR partial pres- 
sure and the wave energy density, respectively, and z' = z — vj:,. Substituting 
PcR in Ea. (j57p and neglecting slow convection with Va in Eq. (j58|) . we arrive 
at the following diffusion equation for / 

dl d KB dl dPcm 

- n -Va- 



dt dz I dz d 



z 



The equation is supplemented with the boundary condition / — > I ism, for 
\z\ — > oo. Outside of the region where Pcr 7^ 0, the last term on the r.h.s. may 
be neglected. The second term may be eliminated by replacing Jexp {Ft) — ^ 
/, /„ exp [rt) dt — > t. However, if P is taken in the form of Ea.ljM]). it is fairly 
small due to the factor yJvg/L ^ 1 . We may simply neg lect it. The solution 
for / and Pcb. (z, t) may be found in an implicit form fsee lMalkov et a l. 2012^ 
for details). However, there exists a very accurate convenient interpolation 
formula that can be represented as follows 



^CR = 



2KB jp) 

3/2 rj— 7 



-3/5 



^5/3 + Pj^^)5/6 g-CV4AsM (65) 



where Lc is the size of the initial CR cloud, C = zj sJvaLc^-, and I?nl 
C (n) DisM cxp (— il), with n being a normalized integrated pressure 



J = — / PcRdz 

KB J 



and DisM is the normalized background difFusivity 

-t^iSM - —j-hsM 

Va^c 

while C ~ 1, for 77 > 1 and C - 77"^ for i7< 1. 

The representation of the solution given in Eq. (j65p is convenient in that 
the function VtPcR (C) does not depend on t, so that the solution can be 
shown for all i, z with only one curve. Figure 1^1 To summarize these results, 
the self-regulated normalized {Vcr. = VuLcPcr/kb (p)) CR partial pressure 
profile Vcr comprises the following three zones (77 ^ 1): (i) a quasi-plateau 
(core) at small z/\/i < ^/TTml of a height ^ (D^i^t)^ ' , which is elevated by 
a factor ^ 77~^ exp (77/2) ^ 1, compared to the test particle solution because 
of the strong quasi-linear suppression of the CR diffusion coefficient with 
respect to its background (test particle) value Disu- Dnl ~ ^iSMexp(— 77) 
(ii) next to the core, where V-Dnl < z/^/i < V-Dism, the profile is scale 
invariant, Vcr ~ 2/z. The CR distribution in this "pedestal" region is fully 
self- regulated, independent of 77 and 7?ism for 77 :^ 1, (iii) the tail of the 
distribution at z/^/t > \/Djsm is similar in shape to the test particle solution 
in ID but it saturates with 77 ^ 1, so that the CR partial pressure is ex 
(Disut)^ exp (— 2;^/47?isMi) , independent of the strength of the CR source 
77, in contrast to the test-particle result which scales as ex 77. Because of the 
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CR diffusivity reduction, the CR cloud half-life is increased and the cloud 
width is decreased, compared to the test particle solution. 

Depending on the functions U (p) and Dism (p) , the resulting CR spec- 
trum generally develops a spectral break for the fixed values of z and t such 
that z'^/t - Dnl (p) ^ Asm cxp (-77). 



9 Summary 

Cosmic rays, being a highly non-equilibrium component, often comprise an 
energy density that is comparable to the ram pressure of energetic plasma 
flows and magnetic fields in astrophysical sources with high energy release 
such as supernova remnants, fast stellar winds, and astrophysical jets of dif- 
ferent scales. CRs may also play a role in the global dynamics of interstellar 
gas in galaxies, in particular, they may support galactic winds. In the pres- 
ence of gravitation, the buoyancy of CRs and magnetic field at galactic scales 
may result in the magnetic Parker instability ( Parkeiiri966l . ll967L IShTjll974L 



IRvu et aDl2003llHanasz et al.ll2009[ ). The local CR diffusion is an important 
factor for the Parker instability to occur. 

The microphysical instabilities discussed above lay the groundwork for de- 
tailed simulations of the global interstellar matter dynamics. In this review we 
addressed the recent progress in understanding of the CR-driven instabilities 
with special attention to non-relativistic shocks. We started with a quasi- 
linear analysis of the growth rates of the instabilities driven by anisotropic 
and inhomogeneous CR distributions. Time dependent nonlinear simulations 
are needed to draw conclusions about the saturation level and the spectra of 
magnetic fluctuations produced by the non-equilibrium CR distributions. We 
used numerical simulations to illustrate the nonlinear dynamics of magnetic 
fluctuations. The CR-driven instabilities are shown to be crucial for model- 
ing particle acceleration sources and the CR escape from the sources into the 
interstellar matter. 
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10 Appendix A 

The dispersion equation Eq. (|20p can be expressed in the elementary functions 
by evaluating Eq. (I22|) and Eq.(l23l) as 
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